EFFECTS OF PERMEABILITY AND VISCOSITY IN LINEAR POLYMERIC GELS 

B. CHABAUD* AND M. C. CALDERERt 

Abstract. Wc propose and analyze a mathcmatieal model of the mechanics of gels, consisting of the laws of balance of 
mass and linear momentum. We consider a gel to be an immiscible and incompressible mixture of a nonlinearly elastic polymer 
and a fluid. The problems that we study are motivated by predictions of the life cycle of body-implantable medical devices. 
Scaling arguments suggest neglecting inertia terms, and therefore, we consider the quasi-static approximation to the dynamics. 
We focus on the linearized system about relevant equilibrium solutions, and derive sufficient conditions for the solvability of the 
time dependent problems. These turn out to be conditions that guarantee local stability of the equilibrium solutions. The fact 
that some equilibrium solutions of interest are not stress free brings additional challenges to the analysis, and, in particular, 
to the derivation of the energy law of the systems. It also singles out the special role of the rotations in the analysis. From 
the point of view of applications, we point out that the conditions that guarantee stability of solutions also provide criteria to 
select material parameters for devices. The boundary conditions that we consider are of two types, first displacement-traction 
conditions for the governing equation of the polymer component, and secondly permeability conditions for the fluid equation. 
We present a rigorous study of these conditions in terms of balance laws of the fluid across the interface between the gel and 
its environment |20| . and use it to justify heuristic permeability formulations found in the literature [39) . [14) . We also consider 
the cases of viscous and inviscid solvent, assume Newtonian dissipation for the polymer component. We establish existence of 
weak solutions for the different boundary permeability conditions and viscosity assumptions. We present two-dimensional, finite 
element numerical simulations to study pressure concentration on edges, in connection with the debonding phenomenon between 
the gel and the boundary substrate upon reaching a critical pressure. 
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1. Introduction. This article addresses mechanical modeling, stability of equilibrium states and analysis 
of boundary value problems of quasi-static gel dynamics. We assume that a gel is an incompressible and 
immiscible mixture of polymer and solvent, and study the coupled system of equations of balance of mass and 
linear momentum of the components. Boundary conditions are of traction-displacement type together with 
statements of the permeability of the gel boundary to the environmental fluid. This work is motivated by 
problems arising in the prediction of the life cycle of body-implantable medical devices. 

Gels consist of crosslinked or entangled polymeric networks holding fluid. In its swollen state, the polymer 
confines the solvent and, in turn, the solvent prevents the gel from collapsing into dry polymer. Gels are 
abundantly present in nature and occur when materials are placed in a fluid environment. Devices such as 
pacemakers, bone replacement units and artificial skin turn into gel when implanted in the body. The materials 
that constitute a device differ in swelling ratio, this causing a build up of stress at the interfaces, as well as at 
the contact between the device and its boundary support. High stresses, above the manufacturer's guaranteed 
threshold, may cause debonding instability leading to device failure. 

The equations that we analyze encode relevant properties of gel behavior such as solvent diffusion, trans- 
port of polymer and solvent, friction and viscosity, elasticity, and time relaxation. They consist of equations 
of balance of mass and linear momentum for the polymer and fluid components, together with the saturation 
{incompressibility) constraint. Assuming that the polymer is an isotropic, elastic solid, the total energy of the 
gel is the sum of the elastic stored energy function of the polymer and the Flory-Huggins energy of mixing. 
The variable fields of the equilibrium problem consist of the volume fractions of the gel components and the 
deformation gradient tensor of the polymer. These fields are not all independent, but they are related by the 
equation of balance of mass of the polymer and the saturation condition on the gel. The problem of con- 
strained energy minimization studied by Micek, Rognes and Calderer in |26j established sufficient conditions 
on the energy and on the imposed boundary conditions that guarantee existence of a global energy minimizer. 
The authors also developed, analyzed and numerically implemented a mixed finite element discretization of 
the equilibrium linear operator. A challenge of the analysis is the presence of residual stress in the reference 
configuration of the polymer. The work yielded numerical evaluations of shear stress at the interface between 
two gels representing bone tissue and the artificial implant. However, gel behavior is inherently a time evolu- 
tion problem due to the combined effects of transport, diffusion and dissipation. This serves as a motivation 
to the work presented in this article. 
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Scaling arguments for gels consisting of polymer melts justify neglecting the inertia terms and analyzing 
the quasi-static system. Specifically, the size of polymer dissipation effects with respect to inertia results in 
the latter being dominant at time scales on the order of 10~^ seconds [Hj- Of course, such time scales are 
negligible for biomedical devices with typical life-cycle of 20 years. From a different perspective, the effect of 
the inertia terms was analyzed by Zhang and Calderer [9]. They studied the free boundary problem of gel 
swelling, in one space dimension, neglecting the Newtonian dissipation of the both gel components. In such 
a case, the governing equation turns out to be a frictional, weakly dissipative, hyperbolic partial differential 
equation [T^], [13]. The scaling that justifies retaining the inertia effects is consistent with the dynamics of 
polysaccharide gel networks found in many living systems, such as in gliding myxobacteria [22j, [30] . 

The model proposed here shares analogies with deformable porous media flow models but has the additional 
feature of accounting for interaction between fluid and polymer through the Flory-Huggins energy. From a 
different perspective, especially challenging multi-component mixtures are used in geology and in oil and 
natural gas recovery models [1], ^28^. In these models, a relevant role is played by the Terzaghi's stress, 
which is the pressure exerted by the fluid in the pores against the stress applied to the rock. Its analog in the 
case of the gel is the pressure in the solvent accounted by the Flory-Huggins energy. The current analysis is 
not immediately extendable to the triphasic models when one of the components is compressible. 

The governing system that we study is motivated by the stress-diffusion coupling model developed by 
Doi and Yamaue |39|, [Mj, [40], [41], and also by Hong et al. :37]. Feng and He [15] studied purely inviscid 
gels with impermeable boundary as governed by the stress-diffusion coupling model. Our model treats the 
polymer component of the gel as a nonlinear elastic solid, and includes, both, fluid and polymer dissipation, 
diffusion, the Flory-Huggins interaction, and accounts for different permeability properties of the interface 
between the gel and the surrounding fluid as well as traction-displacement conditions imposed on the polymer 
boundary. To our knowledge, all these combined effects have not been accounted for in the previous works. 
We develop an existence theory for the governing system of partial differential equations, linearized about an 
equilibrium state. Some of the tools presented in [15 have been applied to our analysis, in the case that the 
solvent is inviscid. One novelty of this work is the derivation of permeability conditions from balance balance 
laws at the interface between the gel and its surrounding fluid, following a model developed for the treatment 
of polyelectrolyte gels [20] . 

We assume that the polymer component of the gel is an isotropic elastic material with Newtonian dissi- 
pation. The equations possess bulk, stress free, equilibrium solutions that are pure expansion or compression. 
We point out that, as for isotropic elasticity, if the energy is convex with respect to the deformation gradi- 
ent, these are the only stress free critical points. However, if the energy is nonconvex, the system may also 
admit non-spherical equilibrium deformations. These states may be consistent with experimentally observed 
pattern structures in polyelectrolyte gels [5], [5] and [IS]- Equilibrium gel states with nonzero stress are also 
relevant in many applications. In particular, reference configurations with residual stress may be counted as 
a special case of the former. Our analysis, addresses the two types of linearization of the governing system, 
first, about stress free spherical deformations, and secondly about equilibrium solutions that satisfy mixed 
traction-displacement boundary conditions. Within this perspective, we may consider the process of device 
implantation as subjecting an originally stress free body to displacement initial and boundary conditions, as 
well as to the environmental stress of the surrounding tissue, and the permeability effects of such a contact. 
The device will no longer be at equilibrium under the newly imposed initial and boundary conditions, and 
a dynamical process will begin at implantation. The second type of linearization is relevant to the iterative 
process of solution of a nonlinear problem. 

We derive restrictions on the constitutive equations that ensure the coercivity of the static operators. 
These are also known as the Coleman and Noll conditions and guarantee the classical stability requirement 
that the stress work be non- negative in every strain ([3ii,, sections 52 and 83). Moreover, we find that the 
procedure of deriving the energy relation brings out the special role of the rotations in the case that residual 
stresses are present. 

We assume that the reference configuration of the gel is that of the polymer network previous to the gel 
formation. Accordingly, the boundary of the current domain is that of the polymer, and it evolves with its 
velocity. The Eulerian formulation of the governing system of the gel and the natural Lagrangian setting of 
solid elasticity of the polymer present challenges to the analyses. These manifest themselves in the derivation of 
the energy relations satisfied by weak solutions of the governing systems. We consider the cases of impermeable 
and fully permeable boundary between the gel and the environmental fluid. The case of a semipermeable gel 
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boundary follows from the former, with some elementary modifications, and so, we omit its presentation. 
Finally, we point out that the conditions for local equilibria employed in the solvability of the time dependent 
problems are also sufficient to guarantee regularity of the weak solutions. 

In related work, we developed and analyzed a numerical method based on finite elements to simulate 
solutions of the models presented here, in two dimensional domains, in the case that the fluid is inviscid. We 
consider a gel sample in the unit square, subject to zero displacement in two opposite edges and to a flxed 
pressure of lO^Pa in the other two. We calculate the stress components under the following criteria: the relative 
scaling of the Flory-Huggins energy with respect to the elastic one, the degree of stiffness, expansion and 
compressibility of the polymer, and the type of boundary permeability. The elasticity modulus of the polymer 
is set at IGPa, consistent with values encountered in device materials. We find that pressure concentrates on 
the edges where the the displacement is held to zero, its values increasing with decreased compressibility, and 
large stiffness. Permeability also promotes stress concentration, but with the interior stress being lower than 
that in the impermeable case. If the pressure at an edge overcomes the debonding threshold, then it would 
detach from its support. The experimental literature reports on values of the debonding pressure for different 
materials ranging from 0.5 to 10 times the elastic modulus /i^, when the value of the latter is of the order of 



The paper is organized as follows. In section 2, we present the balance laws of the gel, the constitutive 
equations and discuss the equilibrium states. Section 3 is devoted to the linearization of the governing equa- 
tions, formulation of boundary conditions, and the study of the local stability of the equilibrium solutions. In 
section 4, we study existence of weak solutions in the case that the fluid is inviscid, and section 5 is devoted 
to the case of viscous solvent. In both sections, the linearization is carried out about uniform dilations or 
compressions. In section 6, we derive the energy law in the case that the governing equations are linearized 
about an arbitrary equilibrium solution. This is the main ingredient in extending the analysis of sections 4 
and 5 to the more general case, and for which we omit the details. The numerical simulations are presented 
in section 7. Finally, in section 8, we draw some conclusions. This work is based on the Ph.D dissertation by 
Brandon Chabaud (10) . 

2. Modeling of Gel Mechanics. We assume that a gel is a saturated, incompressible and immiscible 
mixture of elastic solid and fluid. In the reference conflguration, the polymer occupies a domain fl C R^. The 
solid undergoes a deformation according to the one-to-one, differentiable map 



We let Qt — y{^,t) denote the domain occupied by the gel at time t > 0, and denote F — VxY- We label 
the polymer and fluid components with indices 1 and 2, respectively. A point y G is occupied by, both, 
solid and fluid at volume fractions — (j)i{y,t) and 02 ~ 4'2{y,t), respectively. We let vi = vi(y,t) and 
V2 = V2(y,t) denote the corresponding velocity fields. 

An immiscible mixture is such that the constitutive equations depend explicitly on the volume fractions 
4>i,i = 1,2. We let pi denote the mass density of the ith component (per unit volume of gel). It is related 
to the intrinsic density, 7^, by the equation pi = ^i4n, i = 1,2. Moreover 7^ = constant, i = 1,2 define an 
incompressible mixture. Throughout this section, and unless otherwise specified, the V notation refers to 
derivative with respect to the Eulerian space variable y. The assumption of saturation of the mixture, that 
is, that no species other than polymer and fluid are present, is expressed by the equation 



(In some terminologies, this condition is also known as incompressibility). The governing equations in the 
Eulerian form consist of the balance of mass and linear momentum of each component as well as the chain-rule 
relating the time derivative of the gradient of deformation with the velocity gradient: 



lO^Pa 23 . 



y — y(x, i), such that det(Vxy) > 0, x e 17. 



(2.1) 



4ii + 4>2 = 1- 



(2.2) 



dt 



^ + V • (0,v,) 



0, 



(2.3) 




(2.4) 



dF 

Ik 



+ (vi • V)F: 



(Vvi)^^, 



(2.5) 
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rj > constant, y G fit, i — 1,2. Here Ti is the Cauchy stress tensor of the ith component. The Lagrangian 
form of the equation of balance of mass of the polymer component is 



0i(y(x,i),t) detF(x,t) = M^), 



(2.6) 



where <(/),< 1 represents the prescribed, differentiable, reference volume fraction. Adding up both equations 



in (2.31, and taking (2.2) into account gives 



V • ((/)lVi + (?!)2V2) = 0. 



(2.7) 



The governing system consists of equations ( 2.2 )-( 2.5 1 supplemented with constitutive equations for 71 together 



with prescribed boundary and initial conditions. Rather than analyzing this system directly, we will instead 



consider the set of equations (2.2|, (2.4), (2.5), (2.6) and (2.7|. 



2.1. Gel environment. Upon body implantation, the device becomes immersed in tissue occupying the 
domain Bt, with ilt C Bt- Wc let Tb denote the stress in the surrounding fluid. In f20^, we prescribe governing 
equations for the fluid in Bt as well as balance laws at the interface dilt H dBt- In this article, we adopt the 
simplified assumption that the outside fluid exerts a prescribed pressure % = —PqI, and do not postulate 
balance laws in Bt- Letting v;, denote the velocity field of the outside fluid, we assume that the following 
relations hold on dil.t PUI : 



02(v2 - vi) • n =(vf, - vi) • n := w, 
(vfc - vi)|| =(v2 - vi)|| := q, 

(71 + 7^)n + [p]n =-72^^(1 - -^)n. 

02 



(2.8) 
(2.9) 

(2.10) 



Equation (2.8) is the statement of balance of fluid mass across dflt and (2.10) states the balance of linear 



momentum of the fluid across the interface. Here [p] := Pq — p, where p is the gel pressure at the interface 



limit. It is easy to check that the right hand side of equation (2.10) is the change in linear momentum density 



of fluid crossing a unit area of the interface, and the left hand side represents the total force per unit area 
acting on the fluid. 

Also, following ^20j, we assume that the interface has an intrinsic viscosity, with coefficients 77_l > 
and 77|| > 0, respectively, affecting the fluid crossing it in the normal direction, or moving tangentially to it. 
Specifically, we assume that 



■■=lMi^^' - - b] - n • ( J)n = yy^. 

2 02 02 



Hii 



=(rin) 



m[^2 - vij 



(2.11) 
(2.12) 



2.2. Boundary conditions. We assume that the gel is surrounded by its onw fluid. The boundary 



conditions at the interface Silt consist of the set of equations (2.8|-(2.12|. In the case that the mass inertia of 



the fluid is neglected, they yield two types of boundary conditions, traction on the gel and equations expressing 



the degree of permeability of the gel boundary to its surrounding fluid. First of all, from (2.8) and (2.9), we 
obtain the velocity V{, of the fluid outside the gel but near the boundary: 



Vb • n = (0ivi + 02V2) ■ n, Vfc • q = V2 • q, ondflt 



for any q 7^ such that q • n = 0. Equation (2.10), yields balance of force at the gel-fluid interface 



Pon=iTi+T2-pI)n. 



(2.13) 



(2.14) 



equation (2.8), the flrst one becomes 



Equations ( 2.11 )-( 2.12) are statements of semipermeability of the gel interface. Neglecting inertia, and using 

(2.15) 



,T2 



[p] - n • ( — )n = 77_L02(v2 - vi) • n. 

02 
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Now, substituting equation (2.14) into (2.12), the latter becomes 



7^n||| = 77||(v2 - vi)|| 



(2.16) 



So, the boundary conditions on 9rit consist of equations (2.14), (2.15) and (2.16). 



We now take limits in equations (2.15) and (2.16) as 77||,77_l 

and 



2(v2 - vi) • n = 0, 

,T2 



oo, 0, giving, 

(V2 - Vi)|| = 0, 



[p] - n • ( — )n = 0, and 72n||| = 0, 



(2.17) 
(2.18) 



which correspond to the case of impermeable and fully permeable boundary, respectively. Rather than im- 
posing the traction condition (2.14) on the whole interface d^lt, we will consider mixed traction-displacement 
boundary conditions. That is. 



- Pon={Ti+T2-pI)n, onr, 
u =U, on Fq, 



(2.19) 
(2.20) 



where Fq U F = dU.t, Fq n F = 0, and u = y — yo denotes the displacement vector, with yg representing an 
equilibrium deformation field to be chosen later, and U a prescribed boundary displacement. 

Remark. The condition of semipermeability states the continuity of the force acting on the fluid across the 
interface. A related expression, with 02 = 1 in equation (2.15), is usually found in the literature. 

2.3. Energy, dissipation and constitutive equations. The total energy of system consisting of the 
gel immersed in the environmental fluid is 



1=1,2 

V':=</)iW^(i^) + G(0i,02) 



(2.21) 
(2.22) 



where ■0 denotes the free energy per unit deformed volume of the gel, and W = W{F) and G — G(0i,02) 
represent the elastic energy density of the polymer and the Flory-Huggins energy of the gel, respectively, with 



G(0i, 02) =a0i In 01 -I- 602 ln(?!)2 -I- C0102, 
_ KO KO _ K9 



(2.23) 
(2.24) 



Here K9 denotes the macroscopic energy unit, Vm > is the volume occupied by one monomer, Ni and N2 
represent the number of lattice sites occupied by the polymer and solvent, respectively, and the number 
of monomers between entanglement points; > denotes the absolute temperature, and x > represents the 
Flory interaction parameter [16] and [17]. We let 72' and 7^",« = 1,2 denote the reversible and the viscous 
components of the gel stress, respectively. So, the total stress of the i-component is Ti — + T^" , i — 1,2, 
with T^J — 0. Moreover, we assume that the dissipative stress of each component is Newtonian, 



D(v) := i(Vv + Vv^), 



(2.25) 



i — 1,2; rji > and fii > denote the shear and bulk viscosity coefRcients, respectively, of the components. 
The total stress T of the gel is 



T = Ti+T2. 



(2.26) 



From now on, we will treat (2.7) as a constraint, and let p denote the corresponding Lagrange multiplier. 
Letting e(x, t) denote the internal energy density of the gel and C(x, t) — |(e(x, i) — ip{x,t)) the entropy 
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Fig. 2.1. These are plots of the osmotic pressure of the gel with respect to the polymer volume fraction <p. We point out 
the change of monotonicity of the graphs with respect to x- Such phenomenon corresponds to phase separation in the gel. 



density, the Clausius-Duhcm inequality states that C > holds, for all admissible processes of the gel. This 
allows us to establish the following proposition, which proof is presented in [8J . 

Proposition 2.1. Suppose that the Clausius-Duhem inequality holds for all admissible processes, and let 



{vi,(^i,p} smooth solutions of equations (2.2), {'2.4), [2.6), (2.1), (2.11), (2.12) and (2.25). Suppose that the 



boundary conditions (2.S)-(2.12) are satisfied. Then the following relations hold: 



7T 



a — ttJ, a = (i 

dG _ ac, 



^ dF 



(2.27) 
(2.28) 



Moreover, the dissipation inequality 



dT 
~dt 



ryJD(vi) +^i(V-Vi) + ?7 vi - V2 My 



{lUw"" +v\\H^)dS, 



(2.29) 



holds, where n denotes the unit outward normal to the boundary. 

As a result of the required material frame-indifference, there exists a function W defined on the space of 

symmetric, positive definite tensors, such that W{F) — W{C), where C = F^F. We let V — 2(f) j ^'^'^^ denote 
the second Piola-Kirchhoff stress tensor. We further assume that the polymer is an isotropic elastic material, 
so there is a scalar function w of the principal invariants I — {Ii, 12,13] of C, such that W{C) ~ w{I). In 
this case, the following representation holds (PT*, page 279): 



dW _ 



ai(I)/ + a2{T)C + ao{I)C~\ 
dw dw ^ dw 



dh 



ai 



dh 



1 1 



0L2 



dw 
dh' 



(2.30) 
(2.31) 



The Cauchy stress tensor a = -^F'P{C)F'^ in (2.27) has the form. 



a = 0:(/3o(2:)/ + Pi{I)B + /32(I)S-i), 

2 dw 



^'-Th^'dh 



//, dh 



p2 = -2^/: 



dw 
dh' 



(2.32) 
(2.33) 
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Example 2.3.1. Let us consider a Hadamard material (\1 If . lESjj), 



w{Ii,l2,h) 
a — 'tltll (i/Q _ i^A 



Pi = -y^^, =0, ^0 = J=- 

Vi3 VJc 



03 



Mb'/'/ 



(2.34) 
(2.35) 
(2.36) 
(2.37) 



where a > 0, ai > 0, 03 > 0, g, s, r > 1 are constant. The parameter represents the crosslink density of the 
network. 



Stress free bulk equilibrium states {(j)o, Fo,po) are constant fields satisfying (2.6|, and 

T{{Fo, M - <PoPoI - 0, -(1 - Mpo = 0, 



with T{ as in (2.271. It is easy to check that the former reduce to 

Pif + fiif-^ = ct>'J^A^) - /3o, where 

B = fl, 4^^^jf-'l, />o. 

For the energy in (2.34), equations ( |2.38 1 become vB = kI, with v and k as in (2.361. Equivalently, 



7r(</)) 



+ ailr~asi^)'^ = a^r''{^)^^ 



(2.38) 



(2.39) 
(2.40) 



(2.41) 



with 7r((/)) as in (2.23) and (2.28). We summarize the previous statements in the following: 

Proposition 2.2. Suppose that the gel is isotropic. Then the stress free bulk states are uniform dilations 
or compressions satisfying equation ( 2.3£^ and (2.40). In the case that the energy is given by {2.34), there is 
a unique equilibrium state provided the material parameters satisfy 



- > ai3"-V/ + 030?"^ + — (c + 6 - a) 



2q 



1 



(2.42) 



Moreover, k(0o) > 0, and 4>o > (f)* > 0, where (f>* satisfies K((j)*) = 0. 

Remarks. We point out that the assumptions on the coefficients of w and G are sufficient to guarantee the 
existence of a global minimizer of the total energy under appropriately prescribed boundary conditions [2], 
including those of displacement-traction type. Inequality (2.42) gives insights on material parameter ranges 
that guarantee existence of unique stress free equilibrium states. 

1. Phase separation may occur for material parameters such that k{4>) is non-monotonic. A necessary 
condition for the latter to occur is that c > in (2.24) be sufficiently large, and therefore it corresponds 
to the case that 7r((/)) has negative intervals; in device applications, usually c = .5. This is illustrated 
in Figure 2.1. 

2. Holding c fixed, small values of a > or > may also lead to phase separation. The latter 
correspond to prescribing small shear and bulk moduli. Moreover, the equilibrium value < 0o < 1 
also increases with respect to fj.^ 

We now let 4>a,FQ denote a solution of the boundary value problem 



V • T{{F, (f>) = 0, (pdet F = (pi, p = constant, 
subject to mixed displacement-traction boundary conditions. 



(2.43) 
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3. Linear problems. We now linearize the governing system (2.4), (2.5 1, (2.6) and (2.7) about a partic- 
ular time independent solution {(j)Q, Fq,pq = 0,^1 = = ^2). Relevant special cases include stress-free dilation 
or compression states Fq = /q/, and also non-stress free equilibria. The tilde notation represents perturbations 
from the equilibrium state, which is labeled with 0-super (or sub) indices. Let F := VxU, and 



= (?!)0 -I- 0, F = Fq + F, p = p, Vi = V, V2 = w. 



(3.1) 



The gel domain now corresponds to the reference configuration, 51, of the polymer. First, we calculate the 
linear swelling ratio and polymer volume fraction: 



det F = 00 det ^^o(l + 7(Vxu)) + o{\W^u\^) 
= 0o(l-7(Vx£i)) + o(|Vxiin, 

where 7(Vxu) = tr (Fq^^Vxu). Letting Co = Fq, we denote 



^° =71(00,1-00), 7r° = 



6o,l"0o), i-1,2. 



(3.2) 
(3.3) 



(3.4) 
(3.5) 



The fourth order tensor with components ^ijki corresponds to the elasticity matrix, with the symmetry 
properties 



^ijkl — ^klii — ^iikl — 



^klij 



■jikl 



~ijlk ■ 



(3.6) 



The quantities (3.2)-( |3.5[ ) yield the linearized expressions of the stress tensors. These equations are exact up 
to terms of order o(|VxUp)): 



a =ao + 0o{-i^o7'oJ^o^7(Vx£i) + Foe:(Fo^VxU + Vx£i^Fo)Fo^ 
+ (FoPoVxOi^ + VxOiPoFo^)}, 

7r-7r"-0o(7r?-^°)7(Vxii), 



r =T{ = V 



r,0 



/.o(7r?-^2")/-a(Fo,0o))7(Vxfi) 



+0o{Foe:(Fo^VxU + VxU^Fo)Fo^ + (FoT'oVxU 
The linearized system of equations is 



.~.T 



V^uVoF^)} 



V • (00V -f- (1 - 0o)w) = 0, 

V • Ti - /3(v - W) - 00 Vp = 0, 

V • + /3(v - w) - (1 - 0o)Vp ^ 0, 
ri(Vxii, Vv) = T{iV^u) + ryiD(v) + /ii(V • v)/, 

T2(VxU, Vw) = 772D(w) + /i2(V • w)/. 
Ut = V, 



(3.7) 



(3.8) 



(3.9) 
(3.10) 
(3.11) 
(3.12) 
(3.13) 
(3.14) 



together with (3.3). We point out that the last equation follows from the linearization of (2.5), neglecting 



uniform translations: 



Vxii = (Vvi)i^o- 



(3.15) 



3.1. Stability of Equilibrium Solutions. The conditions that guarantee the stability of the equilib- 
rium states turn out to be also necessary conditions for the solvability of the time-dependent, quasi-static 
problem, that we study in later sections. In order to established such conditions, we first outline the proper- 
ties of the second order tensor V and that of the fourth order one £. Unlike the case of linearizing about a 
stress free state, here we need to include Tr in the analysis of stability. 
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Notation. With the understanding that the quantities that we study are evaluated at Fo,Co = FqFq, hi 
this section, we suppress the 0-notation in the equihbrium solution, and the tilde symbol in the perturbation 
terms, unless explicitly needed. The fourth order elasticity tensor 



j{C)^ Ai6pg6,, + A2{tY CSpg - Cpg)aj + A^idet C)C-^C-^ 

2 



pq ^i] 



da.. 



+0i25ip5jq - aoC^p^Cgj^, Am = ^ ^y^, m = 1,2,3. 



(3.16) 



In the case that C corresponds to a pure expansion or compression, C — f^I, / > 0, we obtain the following 
representations. 



^ijpq =Hf)Sij6pq + 2fi{f)5ipS.iq, with 



rdan , dan ,^ , (9a„ ,2\ r-i 



n=0 



dh dh 



dh 



Moreover, the total linearized stress tensor (|3.8|) becomes 



V = 2/iD(u) + Atr D(u)/, with 

ji = (f>i{2n + fo^iai + a2fo + aafo^), 



(3.17) 
(3.18) 



(3.19) 
(3.20) 
(3.21) 



Likewise, the analog of the fourth order tensor ( 3.16[ ) that combines the elastic and Flory-Huggins effects is 

^ijpq = Hf)SijSpq + 2fl{f)SipSjq. (3.22) 



Proposition 3.1. Let / > 0. Suppose that 

/i(/)>0, 3A(/) + 2/i(/) > 0, 



hold. Then € in (3.22) is coercive, that is, there exists a constant fxo > such that 



holds, for all A G M^^^ . For the Hadamard energy in (2.34), with k and v as in (2.36), we have 



(3.23) 



(3.24) 



(3.25) 



3.1.1. General equilibrium state. We assume that {Fq, (/>o) is a solution of (2.43), and let be as in 
(3.8), with the elasticity tensor £ given by ( 3.16[ ). We say that an equilibrium solution is locally stable if 



r'^(Fo,</)o)-Vv>0, 



(3.26) 



for all sufficiently smooth fields u, v satisfying (3.15). We now derive sufficient conditions for the stability of 



equilibrium solutions. Let us introduce the following notation. 



dif dhdh 



T^l,2- 



(3.27) 



Proposition 3.2. Suppose that w = w(/i,/3) and that relations (2.30)-{2.33) and {3.15 ) hold. Let 
< 00 (x) < 1 and Fo{x.) G A^"^", x G il, be an equilibrium solution. Then 



T'- • Vv = V • Vv 



"dt 



(3.28) 
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where 



H := e:(VxU^ Fa) ■ (VxU^ FoY + S)(VxU; Fq) + tr'{F^'V^u), 



dw , 



^{W^u;Fo) := ^ir ((Fq-^VxU^)' + VxuCq-IVxU^) + ai(FoVxU^ • ^^uF^' + |Vxu|' 

1-01 

Proof. Starting with 

r' • Vv - (^r° + MFo^iFo^^i^ + Vxii^Fo)Fo^ 

+ (i^o^oVxU^ + W^uVoFj)} + - 'r2)^)7(^o"'Vxu)') • Vv, 



straightforward calculations that apply (3.15) give 
d 



{^(S/^u^Fo) ■ (VxU^Fo)) = e:(VxU^Fo + Fo^Vxu) • (VxU^Fo), 



dt 
~dt 

^(f'r0 2tr^(i^o^^Vxu)) = ^.^^^tr (^^o^'Vxu)/ • (VxUi^o^^), 



= (FoT'oVxU^ + VxuT'o-Fo^) • Vv, 



where V is given in (2.30). □ 



Next, we establish coercivity of the operator Ti in (3.29). For this, let us write 

H '.= Hi + 7^2 J 



Hi := C2tr2(Vu^Fo) + C^tr^iF^^S/u) - ao\\7uFo^\^ + ^tr (VuCg-iVu^), 



"0, 



H2 := ^tr(Fo-^Vu^)^+ai(FoVu^ • Vu^^^-^ + |Vu|^). 



(3.29) 
(3.30) 
(3.31) 
(3.32) 



(3.33) 

(3.34) 
(3.35) 
(3.36) 



(3.37) 
(3.38) 

(3.39) 



Lemma 3.3. Let a,b, c be as in (2.23) and < x < 0.5. Then ttJ 2 > is monotonically increasing. 



Moreover if for each Ii > 0, w{Ii,l3) is convex with respect to /s, then C3 > 0. 

This condition on x is satisfied in gels used in device applications. The monotonicity of vrj 2 for this range 
of X is illustrated in Figure 2.1. 

Proposition 3.4. Suppose that the assumptions of Lemma 3.3 hold. Assume that ao < and that 
C2 > 0. Then 



Ui > |ao||VuFo-i|2 + C2tr^{Vu^Fo) + C'str^iF^'Vu) 



(3.40) 



Next, we study the coercivity of H. Let us consider the polar decomposition F ~ RV, where R denotes the 
rotation tensor, and V — \fC . Let Ai, i = 1, 2, 3 denote the eigenvalues of C . Let us denote 



N :-(VxU^)i?, 

r —h N N h — Q^o + Q^i(^» + ^j) ■ / ■ 

i . flj^j I\ ij 1\ ji^ flij . ^ — , I J. 

We now calculate T-L2 using its representation in terms of the eigenvector basis of C, 



(3.41) 
(3.42) 

(3.43) 
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We combine the first term of the right hand side of H2 with the last one on the right hand side of Hi (3.38) 
(which can also be written in terms of Nij). We also combine the mixed products in Tij with the last term in 
■^2 7 upon application of the Cauchy-Schwartz inequality. We now state 

Theorem 3.5. Let (_Fo,<^o) be an equilibrium solution. Suppose that the assumptions of Proposition 3.4 



hold. Furthermore, we assume that C2 > m {'3.27) and 



ai > - max \hij\ 



Then 



n > 



\ao\ 



\VnF-'\' 



+ C2tr^(yu^Fo) + C^tr^iF-^Vu). 



(3.44) 



(3.45) 



Remarks. 

1. Inequalities (3.23) and (3.24) (for a spherical equilibrium state), and the positivity of C2 and C3 in (3.27) 



(for an arbitrary equilibrium state) correspond to the strong ellipticity of the linear operator. Strong ellipticity 
guarantees regularity of the weak solutions of the linear problem. In the case that Fq = /o J, the assumptions 



of Theorem 3.5 imply inequalities (3.23) to hold 



2. The need to separately account for stretch and rotation in the proof of Theorem 3.5 is a signature feature 
of linear elasticity, when the equilibrium state is not stress free. In particular, the theorem applies to the 
linearization about the reference configuration, even if the residual stress is nonzero. In this case, inequality 



(3.44) is identically satisfied. 



3.2. Initial, boundary-value problems. We formulate the governing equations in terms of homo- 
geneous boundary conditions on the displacement field, which also satisfies the only initial condition to be 
specified in the problem. 



u|f=o uo, uolr,, 



t=o- 



(3.46) 



The latter is a compatibility condition with the boundary data at f = 0. Assume that Fq is of class C™, 
for some given integer m > 1, and U e i/™^^/^(ro). We let U denote the extension of U to fl, so that 
l|U||ffm(n) < C'||U||^™-i/2(ro) f^> P- 68]. From now on, we wiU set m = 2. We also assume Pq G H^/'^{dVl) 



and dn(^C^. Then 3P e H^{n) such that P|an 
Let, 



Po and ||P| 



Hi(n) 



<C\\Pq\\hU2 



(an)- 



u = u-U, Uo=Uo-Ut=o, V2 = V2-V and p = p - 
where V and P are defined as follows: 



P, Ut=Vi, 



V = Ut, P = 

V = 0, P is the extension ofPo, 

for impermeable and fully permeable boundary, respectively. The governing system reduces now to 

V • (</.oUt + (1 - 0o)v2) = 

V • ri(VxU, Vvi) - 0oVp - /3(ut - V2) = fi, 

V • r2(VxU, VV2) - (1 - 0o)Vp + p(ut - V2) = f2, 



with 7i and T2 as in (3.12) and (3.131, respectively, and 



ioUt + (l-0o)V), «= -(l-0o)^ 



fi = -V • ri(D(U), D(Ut)) + 0oVP + f3{\5t ~ V) 
f2 = -V • r2(D(V)) + (1 - 0o)VP - /3(Ui - V) 
G = (P - Po)n - r(D(U), D(UO, D(V))n, ondn, 



H ^h + V ■{- 



)f2 :=-V-H, 



(3.47) 

(3.48) 
(3.49) 

(3.50) 
(3.51) 
(3.52) 

(3.53) 

(3.54) 
(3.55) 
(3.56) 

(3.57) 
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Notation. We suppress the superimposed baron the unknown fields, and write (u, V2,p). 



Without loss of generality, in sections 4 and 5, we consider linearization of the original system about 
uniform expansion or compression, Fq = fol- With the help of the energy law developed in section 6, these 
results can be easily extended to the case of a general equilibrium state. 



4. Inviscid solvent. We prove existence and uniqueness of weak solution in the case that the fluid 



component is inviscid. Setting r]2 = and /X2 = in (3.13 1 and solving it explicitly for V2, yields the governing 
system: 



V2 = vi - (Vp + f2) , 

V • (vi - kVp) = H, 
V-r=fi+f2, 

r = r + 'livvi + vvf ) + /ii(v • vi)j. 



(4.1) 

(4.2) 
(4.3) 
(4.4) 



with as in (3.8). We will analyze two cases that correspond to impermeability and full permeability of the 



boundary, respectively. 



4.1. Impermeable boundary. We assume that the boundary of the gel is impermeable to solvent, so 



that the normal component of the vectorial condition (2.171 holds on dfl. This combined with equation (4.1 ) 
reduces to requiring 



(4.5) 



on the pressure. Moreover, following Feng and He |T5j, we define the variable q = W ■ u, which measures the 
volume change of the solid network of the gel. The system of equations can be reformulated as 



g - V • u, 

- V • (kVp) - H 
V • nW^u, Wut,qt) - Vp = fi + f2, 
r(VxU,Vut,qt) = Atr(D(u))/ + /lD(u) +77iD(ut) 



(4.6) 
(4.7) 
(4.8) 
(4.9) 



with A and fl as in (3.251. The quantities i?, fi, f2 and G are as in (3.53 1-(3.57) with P ^ 0. The initial and 



boundary conditions on solutions of this system are 



u|t=o = uo, q\t=n = V • Uq, 
u|ro = O,uo|ro 0, 
T(VxU, VxVi,p, gf)n|r = G, 



(4.10) 
(4.11) 
(4.12) 



together with (4.5 1. In order to prove existence of weak solution of the governing system, we first derive an 



energy law. For this, we multiply (4.8) by Uj and use equations (4.6) and (4.7) and integrate by parts over fl, 
applying the boundary conditions. 



dt .4 2 



;/i|D(u)|2 + Aq2)dx 



1 



{K\S/p\^+,n\^{ut)\^+fiiq^)dx 



G ut 



(fl 



Uf dx. 



(4.13) 
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We introduce the function spaces and notational conventions: 

W{n) = V e H\n) : u;|r„ = 0}, ^ {u e H\n) : a;|ao = 0}, 

'i>o = {qeL^{0,T-L^{n)): [ fqdxdt = 0}, 

Jo Jn 

Jo Jn 

i-T 



W2 = V e i^(0,r;W(17)) : [ f cV-a;dxdt = 0, VcGR}. 



We write uj G iJg to indicate that every component of the vector function a; is a scalar function in H^{il) 
vanishing on the boundary. 

Definition 4.1. {u,q,p) e W x L'^{Q) x H^{V,) is a weak solution if\/ {u},^p,ip) e W x L'^{VL) x H^{9?), 

ipqdx= / (ysV-udx, (4-14) 
n Jn 

qtTpdx+ / KVp-V-4)dyi= / H^] dyi, (4.15) 
n Jn Jn 

/ {( - p + Ag + ^llqt)I + AD(u) + '7iD(ut)} • V{u}) dx 
Jn 

= /" G-u;d5+ / (fi +£2) -wdx, (4.16) 

/ D(u(0))yD(a;),j dx = /" D(uo).y D(u;)„- dx, /" g(0)¥'dx= ( qofdx. (4.17) 
Jn Jo Jn Jn 

Since no boundary conditions are prescribed on p, an inf-sup condition is needed to establish compactness. 
Lemma 4.2. For any positive T < 00, there exists aQ > such that 

it^/M; ->aom\LHo,T;L^n)), ycj,eL'{0,T;L\n)). (4.18) 

Proof: The proof presented here is due to Sayas ^31^ but is a special case of the general LBB condition [7]. 
Since L^{0, T; ^^(ri)) = $0 © R- (that is, $0 is the orthogonal complement of R under the inner product), 
it is clear that the inequality ( |4.18 1 is equivalent to 



\ fn L{(t>o + c)V -wdxdtl , , 

sup °|| ->aomo\\LHo,T;LHn)) + \c\] V(/.o e $0, Vc e R. (4.19) 

v£L2{0.T-W{n)) l|iJ(,Wj||i2(o^T;L2(o)) 



By [18], (4.19) holds if and only if the following are valid: 



1. There exists an ai > such that 



I ffT fn 00 V • w dx df I 

sup 11°/" II > ai||0o||L2(o,T;L2(n)), V0oe$o, 

v,eL^{o,T;W{ny) \\'-'[^)\\L^{o,T-L^{n)) 



2. There exists an a2 > such that 



I fn fo cV • wdxdtl 

W^^^ ^>«2|c|, VceR, 



3. L2(0,T;W(17)) = Wi +W2. 
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Note that since L^{0,T; H^{n)) C L^{0,T;W{n)), the first item holds if 



sup 

wGL2(0,T;Hi(O)) 



|D(w) 



lL2(o,T;L2(n)) 



> aiWMlL^ia.T-L^ifi)), V(/)o e i^(0,T;$o)- 



The latter is a well-known result shown in [19] . In order to prove the validity of the second item, note that 
for w € L^(0, T; W(fi)), J^cV-wdxdt = Jp cw • mis' dt. Thus the result holds if it is possible to find 
awe L^{0,T;W{n)) satisfying w • n 7^ 0, where n denotes the unit outward normal to T. Assuming 

that r is Lipschitz, we take Fi C F with nonzero measure, and a fixed vector m e R" such that for some 
<5o >0, 

m • n(y) > So 

for a.c. y G fi. Choose any ip G C°°([0,T] x F) with ip > 0, suppf C [0,T] x fi, and Jf^ ip > 0. The 
function ip : [0,T] x F ^ R can be hfted to an element w € L^(0, T; whose trace on [0,T] x F is i^j. 

Take w = wme L2(0,T;W(r2)). Then 



f I •w-ndSdt= I I pm- ndSdt> So [ [ pdSdt 
Jo Jr Jo Jfi Jo Jfi 



> 0. 



Hence the second item holds. Finally, to prove item 3, we must show that for any w £ L^(0, T; W(ri)), there 
exist wi e Wi and W2 G W2 such that w = wi + W2. Note that Wi is equivalent to the set of vectors 
in L^(0, T; W(J7)) with constant divergence. Also, W2 is equivalent to the set of vectors in L^(0, T; W(f7)) 
with normal component on F equal to 0. Since L'^{0,T; H^) C yV2, sele ct W 2 G L^(0, T; iJp (fi)) satisfying 
V • W2 = V • w — /o^/o ^ ■ Set wi = — W2. This implies that (4.19) holds and thus completes the 



T\n\ 

proof of the lemma. □ 

We are now ready to prove the following theorem. 

Theorem 4.3. Assume that the hy potheses of Lemma 3.3 hold. 



Let Fr 



SL. 



h, det Fq denot e 



an equilibrium solution satisfying (2.39). Suppose that A and jl are as in (3.20)-(3.21) and satisfy (3.23). 
Assume that for some T > 0, the prescribed boundary conditions satisfy U G H^{0,T; H's (Tq)) and Pq G 
i^(0, T; L^(F)); let Uq G W(0) denote the prescribed initial displacement. Then there exists a unique weak 
solution {u,q,p) to the initial boundary value problem [^.jp - pTil ) that satisfies 



UG L°°(0,r;iji(f^)), ut G L\0,T-H\n)) 
q G L°-{0, T- L^i^)), qt G L^O, T: L^fl)) 
peL^Q,T;H\n)). 

Proof: First of all, we note that the right hand side terms f^, H and H of the governing equations are given 



by (3.53l-(3.57l with P = and Q = V • U. We apply the Faedo-Galerkin method together with the discrete 
version of the inf-sup condition of Lemma 4.2. For this, we decompose W(f2) = Wi(r2) W2(r2), where 
Wi is the set of all divergence-free vectors in W(r2) and W2 denotes its orthogonal complement under the 
inner product /^^ D(w)ijD(vi)ij, for G W(f2). Wi and W2 are both separable Hilbert spaces, so there 

exist sequences of linearly independent smooth functions {w(^)'''}^j^ and {w'^^^''^}^-^ which are dense in Wi 
and W2, respectively. Moreover, the sequence {w^^^^'^, w^^'''^}^-^ forms a linearly independent dense set in 
W(ri). Define cj)'' = V ■ vir(^)''^ for k = 1, . . . ,00. {^'^jfeLi forms a linearly independent dense set in Since 
^ ^ {'^'^IfcLi forms a linearly independent dense set in \1/ as well. For any integer > 1, define the finite 
dimensional Galcrkin spaces 



Wn = span{w(l)'^w(2).'=}f^^, cl>^ = spanj^'^lf^i, - <f 



We now establish the discrete inf-sup condition. By (4.18) and according to JL8^, it can be shown that for the 



same ao as in (4.18) and for each A^ > 1 



sup 

w"eL2(o,T;Wjv) 



1/0 



N 



a 



|D(w 



lL2(o,T;L2(n)) 



> «o||0^||L=(o,T;L2(n)), V<^^ G L^(0,T;$w) 



(4.20) 
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Next, we set up the finite dimensional approximation of tlie problem. We look for (u^, e W^r x <i>Ar x 



^fjv satisfying the following integral relations, for all (w , , -0 ) e Wjv x <I>^ x vf^v: 
Jn Jn 

[ q^ip'^ dx+ [ K\/p^ ■ V^^ dx = - [ n- VV-^ dx, 
Jn Jn Jn 

Ag" + ^igf )/ + mD(u^) + 7?iD(uf )) • V{u;^) dx 



G-u}'" dS+ / (fi + fa) • dx, 



/ D(u^(0)),,D(w^),;, dx= f D(uo),jD(w^),, dx, 
Jn Jn 



q'^0)(l)'' dx = / qQ<j)'^ dx 



(4.21) 
(4.22) 

(4.23) 
(4.24) 
(4.25) 



This leads to a system of linear ordinary differential equations in time for the coefficients of with 
a complete set of initial conditions. So, there exists a unique triple (u^, q^ ,p^) S Wjv x $jv x '^n satisfying 
the system for all t € [0,r]. As in the continuous case, it can be shown that the discrete system has the 
following energy law: 



/ 

Jn 



i(A!D(u^(r))r 



A(g^)^(r))dx + 



i(,.ivp^i^+^iiD(ur)i^ 



= { / G ■ uf dS + / (fi + fa) • uf dx} dt + 



1 



m|D(u^(0))| 



fii {qt ) ) dxdt 
A(g^)^(0))dx. 



(4.26) 



Usin g (4.26), Korn's inequality, the initial conditions (4.24) and (4.25), and the discrete inf-sup condition 
( |4 20[ ), we find that u^, g^, , and are uniformly bounded in L°°{0,T; H^{n)), L°°{0,T; L^{n)), 
l7{0,T;H^{n)), L^[Q,T;L^(n)) and ^^(o, T; iji(f])), respectively. Since T is finite, we find upon passing to 
subsequences, that 

• 3u e H^{0,T;H^{n)) n L°°{0,T;H^{n)) such that ^ u in H'^{0,T; H^{n)) and ^ u in 

• 3q e H^{0,T;L^{n)) n L°°{0,T;L^{n)) such that q'^ ^ q in H^O,T; L^{n)) and q^ ^ q in 
L~(0,r;L2(f])); 

• 3p e L2,^o,r;7?i(f])) such that ^ p in ^^(o, T; iJi(f])). 

As in (Temam it can be shown that the triple (u, q,p) is a weak solution of the system. Finally, uniqueness 
of the weak solution follows from the energy law and the inf-sup condition. □ 

4.2. Fully permeable boundary. In this section, we consider the case where the boundary of the gel 
is fully permeable to its surrounding inviscid solvent. The governing equations consist of (4.6l-(4.9) together 
with (3.53)-(3.57). The initial and boundary conditions are given by (4.10 )-(4.12 ) and (2.18), which for 
7/2 = 0,1^2 — 0, the latter reduces to 



p = 0, on dQ. 



(4.27) 



The energy law has the same expression as in the impermeable case (4.13). 

Definition 4.4. A triple ( u, q,p ) e W x L'^{Q,) x Hq{Q,) is called a weak solution if for all {-w,(j>,i/j) e 
W X L^(il) X i?o(r2) equations hold. We now state the following theorem. 

Theorem 4.5. Assume t hat th e hypotheses of Lemma 3.3 hold. Let F n — fnl , (j)Q — (t),dct Fo de note an 
equilibrium solution satisfying (2.39). Suppose that A and jx are as in (3.20)-(3.21) and satisfy (3.23). Let U 
and Pq &e as Theorem 4-3, and Uq G W(r2) denote the displacement initial condition. Then there exists a 
unique weak solution {u,q,p) of problem (4.5)-(4A^ which satisfies 

uGL^{o,T;HHn)), ut £ l\o,t- H\n)), p e L''{0,T; H\n)), 

qeL°-{0,T;L^n)), qt E L^O,T; L^n)). 
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Proof: We define the function spaces 



Vi = {w e W : V • w = 0}, V2 = {0} U {w e W \ Vi : w\on 
V3 = {0}U[W\(ViUV2)]. 



0}, 



We point out that Vi, V2, and V3 are separable Hilbert spaces. Therefore, there exist sequences {w(^^''^}^j^ C 
Vi, {w'^^''^}^-^ C V2 and {w^'^)''^}^j^ C V3 of Unearly independent smooth functions which are dense in Vi, 
V2 and V3, respectively. The sequence {w^^'''^, w^^^'*^, w^"^)'*^}^]^ forms a linearly independent dense set in 
W. Define ^ y . ^(2),fc ^{3),k ^ y . ,^(3),fe for fc = 1 , . . . , oo . The sequence {0(2),^^ ^(3),^ j.o^^ fo^ms 

a linearly independent dense set in $. Since C $, and it consists of functions which are zero on dft, the 
sequence {4>^^'^''^}'^^i forms a linearly independent dense set in vj/. For any integer > 1, we define the finite 
dimensional Galerkin spaces 



Wtv = span{w 



k ^{3).kxN 



$^ = span{0(2)■^0(3).'^•}^^^, 



N 



span 



It is easy to check that the discrete energy law (4.261 holds as well. By the theory o f line a r diff erential 
equations, for each N there exists a unique {u^,(f'~/p^) £ Wat x x "^n satisfying (4.21 1-(4.25 1 for all 
t € [0,r]. Integrating in time over [0,T] and applying the initial conditions (4.10), and using the discrete 



u m 



q m 



energy law, we conclude that the sequences , q'^ , , q^ and are uniformly bounded in L°°{0, T; H^{n)), 
L°°{0,T;L^{n)), L^{0,T;H\n)), L^{0,T; L^{n)) and L^{0,T] H'^{n)), respectively. Since T < oo, passing to 
subsequences gives 

• 3u e H^{0,T;H^{n)) n L°°{0,T;H^{n)) such that ^ u in H'^{0,T; H^{n)) and 
L~(0,r;i/i(r!)); 

• 3q £ H^{0,T;L'^{n)) n L°°(0, T; ^^(O)) such that q'^ q in H^{0,T; L^{n)) and 
L-(0,r;L2(f])); 

• 3pe L'^{0,T;H^n)) such that ^ p in L^{0,T; H^{fl)). 

The conclusion that the triple {u,q,p) is a unique weak solution of the system follows as in the case of 
impermeable boundary. 

5. Viscous solvent. In this section, we consider the problem of a gel immersed in a viscous solvent. 
That is, we take the viscosity coefficients 77,; > and fii > 0, i ~ 1,2, in the constitutive equations of the 
stress. In contrast with the case of non-viscous solvent, with scalar permeability conditions, these are now 
vector relations, for impermeable as well as permeable boundary. The governing equations are 



V • (0out + (1 - 0o)v2) - h, 

V • ri(VxU, V^ut) - 00 Vp + /3(v2 - ut) = fi, 

V • r2(VxU, VV2) - (1 - </)o)Vp + /3(Ut - V2) - f2, 

Ti = 2/2D(u) + Atr D(u)7 + 77iD(vi) + ^i(V • vi), 

T2 = ?72D(V2) + Ai2(V • V2), 



(5.1) 
(5.2) 
(5.3) 
(5.4) 
(5.5) 



with ji and A as in equations (3.25), and /i, fi, ¥2 and G (shown below) as in (3.53)-(3.57). As in the case of 
inviscid solvent, U, Ut, V and P are extensions of the boundary data satisfied by the original variables u, 
Vi, V2 and p, and subject to compatibility conditions. Given ro,r C dil, Fq n F = 0, and T > 0, initial and 
boundary conditions are: 



u|t=o = uq, u|r„ = 
(ri+r2)n|r = G. 



O,uo|ro = 0, 



(5.6) 
(5.7) 



As in the previous section, permeability conditions on F need to be prescribed as well. The selection of V in 



( 3.53 )-( 3.57) will be made according to the boundary permeability. 



5.1. Impermeable boundary. We now assume that dfl is impermeable to the solvent. Accordingly, 
we require that the vectorial boundary condition (2.17) hold. Following Ladyzhenskaya ([21], Ch.l, Sec. 2]), 
we assume that the initial displacement Uq = uo(x) in (5.6) is continuously differentiable and such that 



V • uo = 0. 



(5.8) 
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For the sake of compatibility, we define 



V = Ut in O. 



The latter together with (5.9) imply that 

/i = -V • (0oUt + (1 - 0o)V) in n. 



(5.9) 



(5.10) 



The governing system consists of equations (5.1|-(5.8| and (2.17). It satisfies the following energy relation: 
J f i[/i|D(u)p + A(V.u)2]dx+ / [r;i|D(u,)P+m(V.u,)^+r72|D(v2)P 



2(V • V2)2 + /3\ut - V2I = / G-UtdS- [ (Vfi • D{ut) + Vfa • ^'(va) dx. (5.11) 

Jr Jn 

The function space of the problem is 

2IJ = {(w\w2) e i7i(r!) X H\n) : w\w2|r„ = 0; - w^lan - 0; 
V-[0owi + (l-(/)o)w2]=O}. 

Definition 5.1. A weak solution is any (u,V2) G 2D satisfying for all (w^,w^) G 2U the equations 
{[/iD(u) + AV • u/ + 77iD(ut) + ^iV • uj] ■ D{w^) + [772D(v2) + M2V • V2/] • T>{w^) 

2 

+/3{ut - V2) • (w^ - w^)} dx = [ G • d5 - [ (Vfi • D(wi) + Vfa • D(w2)) dx, (5.12) 

[ D(u(0)) •D(wi)dx= / D(uo) •D(wi)dx. (5.13) 

We now state the following theorem. 

Theorem 5.2. Suppose that {(j)Q, Fq), X and fi are as in theorem 4-3. Suppose that the viscosity coefficients 
satisfy Hi,r]i > 0,i = 1,2. Assume that for some finite T > 0, U G H^{0,T; H^{n)), g G L^O,T; L^{T) ) 
for a.e. t G [0,r], Uq G Hj(r2), and satisfy relations (5.9)-(5.8) . Let ti,i — 1,2, H and % be as in (3.54)- 
(3.51) with P — 0. Then there exists a unique weak solution (u, V2,p) to the initial boundary value problem 



(5.1)-(5.8) and (2.11) such that 



u G H\0,T;H\n)), V2 G L^{0,T; H\n)), 
peL^{0,T;H\n)). 

Proof: !2IJ is a separable Hilbert space, so there is a sequence of linearly independent smooth functions 
{(w^''^', w^''^)}^-^ C 2n which is dense in 2U. For any integer > 1, define the finite dimensional space 

span{(wl■^w2■'=)}f^,. 



N 



For any integer iV > 1, we seek (u^,vf') G Wjy satisfying for all (w^'^,w^'^) G Wn the equation 
/ ((AD(u^) + AV ■ u^J + 77iD(uf ) + Ml V ■ uf 7) ■ D(w^'^) + (772D(v2^) + ^aV ■ v^J) • D(w'''^)) 

+ / /3(uf - v^) ■ (w^-^ - w^-^) = f G- w^'^ - / Vfi ■ D(w^'^) + Vf2 ■ D(w2'^). (5.14) 
in ir in 



It is easy to assert that there exists a unique (u", v^) G SUjv for all t G [0,r]. Take w 



,1.N 



W 



2,N 



and 

v^. Integrating in time over [0,T] and using (5.14) and standard inequalities, we obtain the energy 



inequality 



N\\2 



A|D(u^^(T))|^ + A(V.u^^(T))1dx+ / / [7^,\Bini 
n Jo Jn 

+r;2|D(v2^)P + M2(V • V2^)2 + /3|uf - V2^|2] dx < C[||D(uo)||i 

+ l|U||^l(0,T;Hl(O)) + l|V|li2(0,T;_f/l(n))]- 



-/ii(V • u^ fdxdt 



(5.15) 



lL2(o,T;L2(r)) 
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From this inequality and the fact that T > is finite, uniform bounds for in H^{0,T;H^{i})), and in 
1/2(0, T;iJi(rj)) follow. These yield the existence of weak limits u e H^{0,T; H^{n)) such that ^ u in 
H'^{0,T;H\n)), and V2 € L^{0,T; H^{n)) such that vf ^ V2 in L^iO,T; H^in)). As in Theorems 4.3 and 
4.5, they are weak solutions of the system. Uniqueness of weak solutions is a consequence of the linearity of 



the problem and the energy law (5.11). This completes the proof of the theorem. □ 



5.2. Fully permeable boundary. We now assume that dfl is fully permeable to solvent, and require 
the linearized form of the boundary permeability condition ( |2.18 1 hold; 



- (1 - 4>o){p - Po)n + 7;2D(v2) + M2(V • V2)n = 0, ondn, 



(5.16) 



where Pq is the hydrostatic pressure of the solvent surrounding the gel. Assuming that Pq € H'2{dVl), we 
denote P S H^{^) its extension to the interior of the domain. The governing system consists of equations 
(5.1)-(5.5| with forcing terms obtained as in ( 3.53 )-( 3.57) by setting V = and letting P be as previously 



mentioned. The initial and boundary conditions are as in (5.6)-(5.7 



and (5.16) 



is letting V — 



Ut . This gives 



Remark. An alternate choice to taking V = in in ( 3.53 )-( 3.57 

V • ((/joUt + (1 — 0o)V) = 0, and corresponds to a class of solutions with no motion of the center of mass of 
the gel, with only the relative velocity present. 

The system satisfies the energy relation ( |5.11[ ). Setting the space of test functions as 

W = {(w\w2) e iji(17) X H^{^) : w\w2|ro 0; V • [^ow^ + (1 - (j)a)-N^] = 0}, 



weak solutions of the system are defined by relations (5.12) and (5.13). 

We now state the following theorem, which proof is analogous to that of the case of impermeable boundary. 

Theorem 5.3. Suppose that {(/)q, Fq), A and jl are as in theorem 4-3. Suppose that the viscosity coefficients 
satisfy ^i,rji > 0,z = 1,2. Then the governing system (5.1)-(5.5), with forcing terms (3.5!^ -{3.51) and 
satisfying initial and boundary conditions (5.6[ )- |5/7| ) and {5.1 6^ has a unique weak solution (u, V2,p) satisfying 
u e iJi(0,T;i?i(r2)), V2 e L^{Q,T;H^{Vl)) and p ^ L^{0,T; H^Q)). 

6. Linearization about non-spherical equilibria. Let us consider the governing system linearized 
about equilibrium solutions (0o,-Fb) that do not necessarily correspond to dilation or compression states. In 
addition, such states may not be stress free. The next proposition establishes an energy law for such systems. 



Proposition 6.1. Suppose that the assumptions of Proposition 3.5 hold. Let £ and T) be as in 
and !13.31 ), respectively. Then smooth solutions of the system ^3~^-{3.13) and (3.15), and (3.1)-(3.8) satisfy 
the energy relation 



d 
dt 



[7yi|D(uOP + Mi(V • utf + r/2|D(v2)P + M2(V • ^2)^ + P\ut - vap] dx 

/ G-UtdS- I (Vfi • I?(ut) + Vf2 • i:'(v2)) dx. □ 
Jr Jn 



(6.1) 



(6.2) 



Integrating the previous relation with respect to i, and using the coercivity properties of £ and D established 
in Proposition 3.4, estimates for u follow: 

^ ^|Vui^„-ip + C2tr2(Vu^Fo) + C,iT\F{;^Vn) 

^[(^(^iW^u^Fo) ■ (VxU^Fo)^ + S)(VxU; Fo)^ + < 2 tr '(F(f ly^u)] dx 



< 



[7?i|D(u,)P + Mi(V • Utf + ?72|D(V2)P + A*2(V • V2)2 + p\ut - V2p] dxdt 



Jn 

T 



[ [ G-\itdSdt- [ [ (Vfi • ^(ut) + Vf2 • i:'(v2))dxdi. 



(6.3) 
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Fig. 7.1. Stress components Uyy and Oxy for he = IGPa, parameters s = 3,q = 1.5, r = 4 and Flory-Huggins scaling 
parameter 10^ Pa. Top row corresponds to fully permeable boundary and bottom row to impermeable. 



With this estimate, the well-posedness of the weak hnear systems, in the cases of non-viscous as weh as viscous 
solvent, and for aU types of boundary permeabihty conditions fohow. This allows us to extend theorems 4.3 
through 5.3 to the more general case of non-spherical equilibria with possible residual stress. 

7. Numerical Simulations. We present two-dimensional numerical simulations of the gel models pre- 
viously analyzed, in the case of a viscous gel immersed in an inviscid solvent, and for, both, impermeable 
and fully permeable boundary. The goal is to investigate concentrations of stress that may lead to failure 
of the device, if critical thresholds are attained. We developed a fully discrete numerical method based on 
finite elements. All simulations have been performed using the DOLFIN library of the FEniCS project [T1I21). 
The equations are linearized about a stress-free swollen or contracted equilibrium state, which is consistent 
with the gel having residual stress. We assume that this state corresponds to that of the device previous to 
implantation. 

The domain of the gel is the unit square 51 = [0, 1] x [0, 1]. Wc construct a uniform mesh of 2048 triangles, 
each with height h = . We take a uniform partition of the time interval and use the backward Euler 
method to discretize the PDE system in time. 

We carry out the non-dimensionalization of the equations according to the following choices of scales: 
• Stresses are normalized by the pres sure scale /ze, the elastic modulus of the polymer (2.34). 
• 



The Flory-Huggins energy density (2.23) is scaled by the factor We set x = 0.5 in (2.231 

We take the time scale as T = ^ sec, where r/i denotes the viscosity coefficient of the polymer, 
set the length scale to L = 1cm. 



We 



We impose mixed displacement-pressure boundary conditions as explained in section 2.2. We assume the part 
of the boundary F = {y = 0} U {y = 1} is subject to a pressure, Pq, that we take to be consistent with the 
arterial pressure: Pq — 10^ Pa. Zero boundary displacement is imposed throughout Tq — {x — 0} [J {x — 1}. 
A normalized initial displacement Uq = (j^ sin(27ra;)/o(l — /o),2/(l — cos(27rx))/o(l — /q)) is imposed in il, 
where /o denotes an equilibrium expansion or compression. We compute stress components, labeling normal 
stresses as a^x £^nd cr^j,, and letting a^y denote the shear stress. The simulations address the following issues: 

1. The ratio of energy scales, ^gj\, ■ We show simulations for the elastic modulus fis = lO^Pa which 
reflects values used in polymer made devices. The scale of the Flory-Huggins energy is taken between 
1 and 10^^. Wc set rji = 10^ Paxsec, which results in a time scale of 0.1 sec. 

2. The degree of stiffness, expansion and compressibility of the polymer as represented by the energy 
exponents s, q and r, respectively. 

3. Type of permeability of the boundary. 

We summarize the findings of our numerical simulations as follows. 
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Fig. 7.2. Stress components ayy and a^y for fiE = IfPPa, parameters s = l,q = 1.5, r = 1.1 and Flory- Hug gins scaling 
parameter 10^ Pa; fully permeable boundary. 

1. The Uyy component presents stress concentration on the two fixed displacement edges, Fq, for the 
whole range of parameters that we tested. The components Uxy and Uxx (not included here) show 
corner concentration. In the case of impermeable boundary, the stress distributes almost uniformly 
across the domain, showing higher values than in the permeable case. Gels with permeable boundary 
show a low stress profile in the interior of the domain, with stresses concentrating on Fq. 

2. The boundary stress concentrations of the ayy component may trigger debonding upon reaching a 
experimentally determined threshold value [55]. 

3. For a given set of parameters, stresses in the case of impermeable boundary are higher than their 
permeable counterparts. This reflects the fact that, in a gel with fully permeable boundary, exchange 
of solvent takes place across the interface causing some stress relaxation. Moreover, in the ideal case 
of pure permeability, the fluid exchange takes place without loss of energy. 

4. We have performed simulations with values of s ranging from s = 1 (Neo-Hookean material) to s = 3 
(hard rubber), and for values of r ranging from r — 1.5 (high compressibility) to r = 4. We found 
that raising either of these exponents by 1, it may increase the stresses by at least by one order of 
magnitude. 

5. The stress values, as represented by their maximum and minimum absolute values, show a decreasing 
pattern with the increase of the Flory-Huggins energy scaling with respect to the elastic one, reflecting 
softenning of the material. 

6. The stresss distribution shown in the figures correspond to time equal to one hour. Calculations 
done for the same data after one day, show stress values in the same order of magnitude as the ones 
presented here. 

7. Whereas the values of stresses shown in Figure 2 may be near the debonding pressure threshold, 
those in Figure 1 may have already crossed it. The experimental literature reports on values of the 
debonding pressure for different materials and loading conditions ranging from 0.5 to 10 times the 
elastic modulus [23] . 

8. Conclusions. We analyzed a model of the dynamics of gels that addresses inviscid and viscous sol- 
vent and polymer, permeability and traction-displacement boundary conditions, elasticity and diffusion. In 
particular, we focused on the linearized system about relevant equilibrium solutions and derived conditions 
for the solvability of the time dependent problems. These are also conditions that guarantee local stability 
of the equilibrium solutions of stress-free dilation and compression states as well as general equilibria, that 
is, solutions of traction-displacement boundary value problems of nonlinear elasticity. Although we proved 
well-posedness of the solutions of the time dependent equations linearized about dilation and compression 
states, the energy laws that we derived would allow us to extend the results to the more general linearized 
equations in a straight forward manner. In particular, the latter includes reference states with residual stress. 

The analysis developed in this article answers specific questions arising in applications. Indeed, the as- 
sumptions ensuring stability of equilibrium solutions, and the subsequent well-posedness of the time dependent 
problem are formulated in terms of the parameters of the elastic and Flory-Huggins energies and their relative 
scale. Although this is far from sufficient to identify a material for a specific application, it does provide a 
criteria to eliminate materials for which instability would occur. This would have an immediate effect on 
reducing the number of costly and time consuming experiments to test a certain material for application by 
as much as 50 percent |25| . In addition to the stability characterization of the parameters, the numerical sim- 
ulations provide data that indicate whether gel pressure has reached the debonding threshold. Furthermore, 
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the choice of viscosity and drag coefficients determine the time of relaxation of a disturbance. 

The simulations presented in the paper accurately address boundary conditions encountered in device 
applications, as well as values of elastic modulus ^.e = IGPa and Flory-Huggins parameters of realistic 
device materials. However, the domains that we use are two-dimensional and so, cannot represent realistic 
shapes of devices. Another important feature not addressed in the current research is the stress concentration 
phenomenon at the interface between two different materials of the device. Development of numerical tools 
based on Discontinuous Galerkin methods is currently underway to simulate actual devices more accurately. 

The stress corner concentrations that we found have also been observed in gel membrane experiments on 
drug delivery devices. In this case, though, the presence of ions significantly magnifies the effect |32| . 

From a different point of view, a better understanding of the debonding phenomenon is needed, perhaps 
appealing to the problem of cavitation and cavity propagation. Experimental work on debonding also brings 
out the viscoelastic aspects of the phenomenon, so its treatment may require the adoption of viscoelastic stress 
strain laws [55] . 

The analysis presented here can be extended to treating triphasic models developed in the study of drug- 
delivery devices [34], [35]. However, this extension is not straightforward since the laws of balance of mass in 
the latter case are significantly more challenging. 

From the point of view of analysis, one goal of the forthcoming work is to study the nonlinear problem 
within the context of the Oldroyd-B models of nonlinear elasticity. We point out that Sections 3.1.1 and 6 
deal with the linearization of the system about arbitrary equilibrium solutions. This provides a necessary 
ingredient in the time discretization of a nonlinear model. 
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